Solutions of a three-dimensional multi-term fractional anomalous solute transport model for contamination in groundwater

The study presents a meshless computational approach for simulating the three-dimensional multi-term time-fractional mobile-immobile diffusion equation in the Caputo sense. The methodology combines a stable Crank-Nicolson time-integration scheme with the definition of the Caputo derivative to discretize the problem in the temporal direction. The spatial function derivative is approximated using the inverse multiquadric radial basis function. The solution is approximated on a set of scattered or uniform nodes, resulting in a sparse and well-conditioned coefficient matrix. The study highlights the advantages of meshless method, particularly their simplicity of implementation in higher dimensions. To validate the accuracy and efficacy of the proposed method, we performed numerical simulations and compared them with analytical solutions for various test problems. These simulations were carried out on computational domains of both rectangular and non-rectangular shapes. The research highlights the potential of meshless techniques in solving complex diffusion problems and its successful applications in groundwater contamination and other relevant fields.


Introduction
Fractional calculus is a mathematical discipline exploring fractional derivatives and integrals and their properties and applications.By encompassing non-integer orders, it extends conventional calculus, enabling the examination of intricate physical phenomena beyond the scope of integer-order models.This field finds practical utility in engineering and various applied sciences, including finance, electromagnetism, signal processing, control theory and mechanics [1][2][3][4][5][6][7][8].The concept of fractional calculus allows the order of a partial differential equation to vary with respect to time or space.As a result, the application of fractional calculus in differential equations has been broadened significantly.Fractional differential equations (FDEs) represent a broader category of conventional differential equations, accommodating derivatives of real or complex orders.These equations find extensive application in various areas, such as fluid mechanics, heat transfer, and electromagnetism [9][10][11].Fractional partial differential equations (FPDEs), specifically, offer distinct advantages by providing a more accurate depiction of real-world phenomena that conventional DEs cannot adequately capture.Notably, FPDEs have been employed to model complex and nonlinear behavior in non-Newtonian fluids [12,13].Multi-term fractional order systems exhibit numerous advantages over singleterm fractional order systems.Firstly, they enable more accurate modeling of intricate physical processes and phenomena due to their capability to encompass a broader spectrum of behaviors and dynamics [14].Secondly, multi-term fractional order systems are particularly wellsuited for control applications, as they provide improved stability margins and incredible robustness.Consequently, they effectively help disturbances and uncertainties in the system, resulting in improved performance and heightened reliability [15].Moreover, multi-term fractional order systems offer enhanced efficiency and effectiveness in addressing engineering challenges, rendering them a valuable instrument across various applications, such as biomedical engineering, control systems, and signal processing [15][16][17].In this study, we numerically simulate a three-dimensional (3D) multi-term time-fractional mobile-immobile diffusion equation (FMIDE), which is given as [18].
Vð� r; tÞ ¼ g 1 ð� r; tÞ; � r 2 @O; ð3Þ In the context of the given problem, r denotes the Laplacian operator, and μ and α are known constants.Additionally, Fð� r; tÞ and g 1 ð� r; tÞ are specified functions.Furthermore, the Caputo fractional derivative operator @ g k @t g k applies to the function Vð� r; tÞ, defined by Caputo [19] as follows whereas Γ(.) represent the gamma function.The stability and convergence analysis of the fractional immobile/mobile transport model ( 1)-( 3) has been investigated in previous works [18,20].The fractional diffusion equation characterizes anomalous diffusion phenomena, wherein particles do not disperse uniformly over time but instead reveal long-range correlations, nonuniform spreading, and other intricate behavior.It elucidates such phenomena through the extension of the classical diffusion equation to incorporate fractional spatial and temporal derivatives.The mobile-immobile fractional diffusion equation, a particular instance of subdiffusion phenomena, has been extensively employed across various domains for simulating particle transport within heterogeneous media, where particles exhibit varying rates of movement between immobile and mobile regions [21].This equation shows distinct advantages in the characterization of anomalous solute transport phenomena within groundwater and porous media.Specifically, it effectively captures non-uniform transport behavior resulting from the presence of reactive surfaces and discrete flow pathways [22].The movement of dissolved substances in liquids through porous media like groundwater or soil, is mathematically described by a model that accounts for the complexity influenced by diffusion, advection, and chemical reactions.This process is not straightforward or predictable [23,24].
The anomalous solute transport model has proven to be highly versatile and applicable in various domains.These applications include predicting the movement of contaminants in groundwater and devising effective remediation strategies [25].It also aids in understanding water movement in aquifers and surface water systems, optimizing pumping and recharge rates, and preventing aquifer depletion.Furthermore, the model facilitates the prediction of solvents and contaminant movement in soil, which is crucial for identifying effective remediation approaches.Additionally, it contributes to understanding carbon dioxide migration in subsurface formations, assisting in the development of efficient carbon sequestration strategies to mitigate climate change [25].Researchers have explored the utilization of the mobile and immobile technique to address vehicular challenges in both saturated and unsaturated environments [23].Gao et al. [24] presented experimental findings that showcased the efficacy of segregating the liquid phase within a porous medium into two distinct regions: mobile and immobile.This novel approach facilitates the representation of scale-dependent dispersion phenomena in heterogeneous porous media, thereby enabling a more comprehensive depiction of reactive solute transport processes.
Currently, there has been a notable increase in the investigation of FPDEs which are notoriously difficult to solve exactly.Consequently, a substantial portion of research efforts has been directed towards the development of numerical techniques for their solution [26][27][28][29][30]. Several widely employed approaches encompass employing finite difference or finite element methodologies, alongside fractional calculus techniques such as the Riemann-Liouville or Caputo operators.Researchers are currently focused on developing more efficient numerical methods such as variational iteration methods [31], spectral methods [32] and radial basis function (RBF) based meshless methods (MM) [33].
The utilization of meshless techniques has garnered significant attention in the scientific community due to their special capability to handle intricate geometries and yield accurate solutions without necessitating mesh generation.These numerical approaches do not need fixed mesh; instead, they utilize uniform or scattered data points to discretize the computational domain.Notably, RBF based MM have emerged as a prominent sub-type.These methods employ mathematical functions to effectively interpolate data values at scattered points, providing reliable approximations for solving PDEs [34][35][36][37][38][39].RBF-based meshless methods offer a significant benefit in terms of their ease of implementation and adaptability to handle problems involving irregular geometries.Moreover, RBFs exhibit a 'global' characteristic, wherein the shape functions utilized for approximating solutions span the entire computational domain.Despite the usefulness of RBF-based meshless methods in solving diverse PDEs in engineering and science fields, some challenges exist.One such challenge pertains to the selection of suitable values for RBF parameters, specifically shape and size, which significantly impact the accuracy and efficiency of the methods.Moreover, another important limitation is the substantial computational cost associated with tackling large-scale problems, as these approaches involve solving dense systems of linear equations.Nevertheless, RBF based MM persist as a popular tool for PDE solving [40][41][42].
Local RBF methods [43] have been introduced as a viable approach to address the challenges encountered by global RBF methods.The key advantage of local RBFs lies in their ability to exhibit compact support centered around each individual data point.Consequently, the resulting sparse matrices are better conditioned, facilitating the selection of an appropriate shape parameter and leading to a more accurate and efficient solution of the linear equations.As a result, local RBF methods tend to exhibit better accuracy and computational efficiency in large-scale problems.

Motivation
The difficulty in solving fractional partial differential equations analytically motivates researchers to explore efficient numerical alternatives.PDEs play a crucial role in real-world applications such as engineering, physics, and finance, but finding exact analytical solutions is often impractical due to their nonlinearity and complexity.Consequently, various numerical methods have been developed and evaluated, tailored to specific problem characteristics and computational requirements.Traditional finite difference and finite element methods are widely used, but they may face challenges when dealing with irregular geometries, complex domains, and moving boundaries.In contrast, meshless methods present an attractive alternative as they do not require a pre-defined mesh and can efficiently handle complex geometries and unstructured domains.This adaptability makes meshless methods particularly suitable for fluid dynamics, structural mechanics, and data-driven modeling.This article introduces a meshless numerical scheme for solving PDE models using RBFs to compute spatial derivatives, ensuring accurate representation of unknown functions in higher dimensions.Additionally, the temporal direction is discretized using the Caputo derivative definition.The proposed meshless approach offers several advantages over traditional methods, including eliminating the need for a structured grid, ease of implementation for complex domains, and enabling seamless extension to multidimensional problems.Furthermore, the scheme demonstrates high accuracy and numerical stability, which are crucial for reliable and robust simulations.

Fractional calculus: Analyzing the theoretical foundations of a time discrete scheme
In the context of discretizing the time variable, we begin by introducing key preliminary concepts from functional analysis.Fractional derivatives, an essential aspect of fractional calculus, hold a significant role in this process.To better understand their significance, we outline the fundamental fractional operator definitions that are commonly employed.

Definition 3
The Atangana and Baleanu fractional operator [46] ABC a Definition 4 He's fractional operator [47] @ g k Vð� r; tÞ Introduction to applied functional analysis: A preliminary overview In this regard, one can obtain o ; Next, we introduce the elucidation of inner product within a Hilbert space. ðV; The Sobolev space X 1,p (I) is said to be Additionally,within the scope of this manuscript, we introduce the ensuing inner product alongside the corresponding energy norms within the function spaces L 2 and H 1 .
respectively.Let us define J ¼ T M be the mesh size in time, and Lemma 1.Let us suppose η(t)2C 2 [0, T] and 0 < γ k < 1, then it holds that and jR n j � Proof.Sun et al. [48].
Plugging RBF ψ(kr − r p k) in ( 9), we obtain c ðmÞ ðk r h À r p kÞ ¼ l ðmÞ hk cðk r hk À r p kÞ; p ¼ h1; h2; . . .; hn h ; ð10Þ In this case, the following matrix expression of Eq (10) is obtained by employed the inverse multiquadric (IMQ) RBF cðk r hk À r p kÞ ¼ 1= ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi  for each k = i1, h2, . .., hn h .Eq (11) can be written as As reported in [49], it has been established that the matrix A n h is invertible.Utilizing (13), we can derive the following result Eqs ( 9) and ( 14) implies where The following elucidates the methodology concerning the two-dimensional scenario.
r ðmÞ k Vðr hk ; s hk Þ; h ¼ 1; 2; . . .; N 2 ; ð17Þ The coefficients r ðmÞ k and Z ðmÞ k (k = 1, 2, . .., n h ) can be determined through the following process Similarly, the identical procedure can be extended to the three-dimensional scenario.An approximate semi-discretized formulation for the mathematical model described by (1) with associated initial and boundary conditions is developed herein by employing the suggested meshless methodology.
Let D denote the sparse coefficient matrix derived from the recommended MM approximation.The vectors b and k of size N×1 represent the boundary and initial conditions, respectively, of the given problem.

Time-stepping schemes
The time-fractional Caputo derivative for 0 < γ k � 1 is @ g k Vð� r;tÞ @t g k , which can be written as Taking into account M + 1 equidistant time levels τ 0 , τ 1 , . .., τ M within the interval [0, τ], where the time step is denoted as J and defined as t n ¼ nJ for n = 0, 1, 2, . .., M, we propose a first-order finite difference scheme to approximate the time fractional derivative term.
ð23Þ where @Vð� r;z j Þ @z , is approximated as follow @Vð� r; The more precise form is

Formulation of θ-weighted scheme
The methodology for the time discretization of a 1-term time-fractional order utilizing the θweighted scheme along (24) to approximate the model given by (1) (taking μ = 1), we have we get similarly for n = 0 After implementing the proposed meshless technique (discussed in Section), ( 27)-( 28) lead to Let L represent the weight matrix of the differential operator L, and I denote an identity matrix.Eqs ( 29) and ( 30) reduces to Crank-Nicolson scheme for y ¼ 1  2 .Likewise, this process can be iterated for 2-, 3-, and 5-term time-fractional derivatives.

Numerical simulation and discussion
This section presents a comprehensive analysis of the accuracy and applicability of the suggested meshless technique for approximating the numerical solution of the underlying problem given in (1).The proposed method utilizes IMQ with shape parameter value c = 1 and its effectiveness is assessed for different time fractional orders, including 2-term, 3-term, and 5-term equations.Furthermore, the proposed method is subjected to testing on both non-rectangular and rectangular domains.The temporal step size is set to J ¼ 0:0005, and the spatial domain is defined as [0, 1], unless specified otherwise.The accuracy of the suggested method is evaluated utilizing the following criteria: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where V is the exact solution.
Test Problem 1 The closed-form solution for the model (1), with μ = α = 1 is Vð� r; tÞ ¼ t 2 cosð2prÞ cosð2psÞ cosð2pzÞ; � r ¼ ðr; s; zÞ 2 O; The presented numerical results for Problem 1 are obtained using the recommended meshless approach and are shown in Table 1.The approach uses various parameters, including the number of nodes N, temporal step size J, and fractional-orders i.e., γ 1 = γ 2 = 0.5 for 2-term, γ 1 = γ 2 = γ 3 = 0.5 for 3-term, and γ 1 = γ 2 = γ 3 = γ 4 = γ 5 = 0.5 for 5-term.Also, the final time τ is set at 0.5, while the accuracy is evaluated using the MaxE and RMS.The results indicate that the suggested meshless approach provides better accuracy.Moreover, the results in Table 1 reveal that as the number of nodes and the number of terms in the time-fractional orders increase, the accuracy improves.Numerical results for various combinations of τ, J, and γ 1 = In Table 4, we display the results of numerical simulations concerning Problem 2 under different values of N, while fixed value of τ = 0.5 and γ 1 = γ 2 = 0.5 for the 2-term case, γ 1 = γ 2 = γ 3 = 0.5 for the 3-term case, and γ 1 = γ 2 = γ 3 = γ 4 = γ 5 = 0.5 for the 5-term case.Similarly, in Table 5 numerical results for the number of J and τ, while keeping the nodes N = 15 and γ k = 0.5 fixed.The value of γ k is chosen as γ k = γ 1 = γ 2 for the 2-term case, and as γ k = γ 1 = γ 2 = γ 3 and γ k = γ 1 = γ 2 = γ 3 = γ 4 = γ 5 for the 3-term and 5-term cases, respectively.It should be noted that accurate numerical results have been obtained in all these cases.Figs 3 and 4 illustrate a comparative examination of solutions using 2-term, 3-term, and 5-term fractional orders.In Fig 3, we present a comparison between numerical and exact solutions, accompanied by their respective absolute errors.Whereas Fig 4 displays the results obtained through the 5-term fractional order approach.The observed results strongly suggest that the proposed method yields highly accurate results.
The meshless method stands out for its capability to tackle problems in non-uniform geometries.Unlike conventional techniques, it doesn't rely on a predetermined mesh structure, which is especially useful when dealing with intricate and irregular shapes.Moreover, this approach is versatile and adaptable, as it can effectively solve problems with scattered data points without needing any connectivity information.The effectiveness of this feature has been demonstrated on the computational domains [43] illustrated in Figs 5 and 6, with the corresponding results.These results indicate that the proposed meshless method produces precise outcomes, even in challenging scenarios.

Conclusion
The study presents a meshless computational approach for simulating the 3D multi-term time-fractional mobile-immobile diffusion equation in the Caputo sense.The methodology combines a stable Crank-Nicolson time-integration scheme with the definition of the Caputo derivative to discretize the problem in the temporal direction.The spatial function derivative is approximated using the inverse multiquadric RBF.This integration yields a sparse linear system of equations, resulting in significant reductions in computational expenses and execution time.To assess the accuracy of the proposed approach, three distinct error norms were computed.Furthermore, the effectiveness of the method was demonstrated by applying it to three different test problems with irregular computational domains.The obtained results were subjected to rigorous evaluation, and their accuracy and efficiency were illustrated through the presentation of tables and figures.Notably, the versatility of the proposed technique enables its adaptation for various complex fractional partial differential equations with minimal modifications.

γ 2 =
0.5, γ 1 = γ 2 = γ 3 = 0.5, and γ 1 = γ 2 = γ 3 = γ 4 = γ 5 = 0.5 with N = 15 nodes are presented in Tables2 and 3.The results indicate that accuracy has been reasonably improved in this scenario.Figs 1 and 2 display a comparative analysis of 2-term, 3-term, and 5-term solutions.Fig 1 presents a comparative analysis between the numerical solutions and the corresponding exact solutions, along with their respective absolute errors.Additionally, Fig 2 exhibits the results obtained for the 5-term fractional order.The observations drawn from these figures strongly support that the recommended method yields accurate results.Test Problem 2 The closed-form solution for the model (1), with μ = α = 1 is Vð� r; tÞ ¼ expðÀ tÞ sinðprÞ sinðpsÞ sinðpzÞ; � r ¼ ðr; s; zÞ 2 O; Fig 7 displays a comparative analysis of the numerical and exact solutions.The absolute error between these solutions is also given in the figure.The figure indicates that the proposed method provides accurate results.

Fig 4 .
Fig 4. For Problem 2, the results for 5-term fractional derivatives are presented, indicating a comparison between exact and numerical outcomes along with the absolute error.https://doi.org/10.1371/journal.pone.0294348.g004

Fig 7 .
Fig 7.For Test Problem 3, the results for 5-term fractional derivatives are presented, indicating a comparison between exact and numerical outcomes along with the absolute error.https://doi.org/10.1371/journal.pone.0294348.g007 hn h ðr h Þ ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } p ðr k Þ ¼ cðk r k À r p kÞ; p ¼ h1; h2; . . .; hn h ; ð12Þ